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Abstract 

We derive variational equations to analyze the stability of synchronization for coupled near- 
identical oscillators. To study the effect of parameter mismatch on the stability in a general 
fashion, we define master stability equations and associated master stability functions, which are 
independent of the network structure. In particular, we present several examples of coupled near- 
identical Lorenz systems configured in small networks (a ring graph and sequence networks) with a 
fixed parameter mismatch and a large Barabasi- Albert scale-free network with random parameter 
mismatch. We find that several different network architectures permit similar results despite various 
mismatch patterns. 
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I. INTRODUCTION 



The phenomena of synchronization has been found in various aspects of nature and 
science[T3]. Its apphcations have ranged widely from biologyjH [TU] to mathematical 
epidemiology [Hj, and chaotic oscillators [2], to communicational devices in engineering [3], 
etc. With the development of theory and application in complex networks [12J, the study 
of synchronization between a large number of coupled dynamically driven oscillators has 
become a popular and exciting developing topic, see for example [HI [151 [HI UHl [El 120] . 

To model the coupled dynamics on a network (assumed to be unweighted and undirected 
and connected throughout this paper), we consider, for i = 1,2, ...,N: 

N 

Wi = fiwi, fii)- gJ2 LijHiwj) (1) 
i=i 

where Wi G 3?™" is used to represent the dynamical variable on the ith unit; / : 3?™ xW ^ Sft*" 
is the individual dynamics (usually chaotic dynamics for most interesting problems) on i 
and fii G 3?^ is the corresponding parameter; L G 9?^^^ is the graph laplacian defined by 
Lij = — 1 if there is an edge connecting node i and j and the diagonal element La is defined 
to be the total number of edges incident to node i in the network; H : 3?™" 3?™ is a 
uniform coupling function on the net; and (? G 3? is the uniform coupling strength (usually 
> for diffusive coupling). The whole system can be represented compactly with the use of 
Kronecker product: 

w = f{w,fi)-g-L0H{w) (2) 

where w= {wj,wj, ...wJjY' is a column vector of all the dynamic variables, and likewise for 
^ and /; and ® is the usual Kronecker product [1]. 

The majority of the theoretical work has been focused on identical synchronization where 
maxjj ||ti'i(t) — u^j(i)|| ^ as t ^ oo, since it is in this situation the stability analysis can 
be carried forward by using the master stability functions proposed in the seminal work [8] . 
However, realistically it is impossible to find or construct a coupled dynamical system made 
up of exactly identical units, in which case identical synchronization rarely happens, but 
instead, a nearly synchronous state often takes place instead, where maxjj ||wj(t) — ^^(t)!! < 
C for some small constant C > as t cxd. 

It is thus important to analyze how systems such as Eq. ([T| evolve, when parameter 
mismatch appears. In ^17j, similar variational equations were used to study the impact of 
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parameter mismatch on the possible de-synchronization. To study the effect of parameter 
mismatch on the stabihty of synchronization, and more specifically, to find the distance 
bound C in terms of the given parameters in Eq. ([T|, we derive variational equations of 
system such as Eq. Q and extend the master stability function approach to this case, to 
decompose the problem into two parts that depend on the individual dynamics and network 
structure respectively. 



II. THEORY: MASTER STABILITY EQUATIONS AND FUNCTIONS 
A. Derivation of Variational Equations 

When the parameters /ij of individual units in Eq. ([T]) are close to each other, 
centered around their mean /i, the coupled units Wi are found empirically to sat- 
isfy m.a,'Xij\\wi{t) — Wj{t)\\ < C for some C > as t ^ cxd, referred to as near 
synchronization[17\. When such near synchronization state exists, the average trajectory well 
represents the collective behavior of all the units. The average trajectory w = Y^iLi of 
Eq. ([T]) satisfies: 

N N N 

1=1 1=1 j=l 

1 ^ 

i=l 

since J2iLi ^ij = by the definition of L. The variation rji = Wi — w of each individual unit 
is found to satisfy the following variational equation: 

N 

Vi = Dwfiw, fi)Vi - LijDH{w)r]j + D^,f{w, fi)5ni, (4) 

where ft = YliLi l^i and SfXi = Hi — ft; and D.^^j represents the derivative matrix with respect 
to w and likewise for and DH. The above variational equations can be represented in 
Kronecker product form as: 



ri 



In (S) D^f - g- L(g)DH r)+ In® DJ 6fi, (5) 



where T]i are stacked into a column vector 77 and likewise for 
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B. Decomposition of the Variational Equations 



Since we are deahng with undirected graph, the associated L is symmetric and positive 
semi-definite, and thus L is diagonahzable: L = PAP^[T], where A is the diagonal matrix 
whose ith diagonal entry A, is the ith eigenvalue of L (arranged in the order Ai < A2 < ... < 
Xn)', and P is the orthogonal matrix whose ith column Vi = {vi^i, VN,iY' is the normalized 
eigenvector associated with Aj, and all these Vi form an ortho normal basis of 3?^. Note that 
because of Yl!i=i Lij = 0, we always have Ai = with vi = 1)"^; and since we have 

assumed that the graph is connected, the following holds: Ai = < A2 < ... < A^r. 

We may uncouple the variational equation Eq. ^ by making the change of variables 



(6) 



or more explicitly, for each i, 



Ci = Vi^iVi + V2,ir]2 + ... + VN,ir]N, 



(7) 



to yield: 



(8) 



where <^ = ((1 , (Jr)'^ ■ Note that since YliLiVi = Y^i^ii^i — ^) = 0, and Vi = 
^(1, 1)^, the following holds: Ci = 0, by Eq. Q. 

Note that since the transformation <^ = {P^ Im)'n is an orthogonal transformation, 
||<^|| = llryll with the choice of Euclidean norm. In other words, for ||.|| being the usual 
Euclidean distance, we have: 

N N 

Eiioir-Eii^^ir- (9) 

1=1 1=1 

The homogeneous part in Eq. ([s]) has block diagonal structure and we may write for each 
eigenmode {i = 2,3, ...,N): 

N 



(10) 



The vector J2j=i is the weighted average of parameter mismatch vectors, weighted 

by the eigenvector components associated with Aj, and may be thought of as the length of 
projection of the parameter mismatch vector onto the eigenvector Vi. 
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C. Extended Master Stability Equations and Functions 



The variational equation in the new coordinate system Eq. (10) suggests a generic 



approach [H] to study the stabihty of synchronization for a given network coupled dynami- 
cal system investigating on the effect of Aj and (^YljLi'Vj^iSfij^ on the solution of Eq. (10). 



We define an extended master stability equation [22j for near identical coupled dynamical 
systems: 



e= D^f-a-DH +DJ-ij 



(11) 



where we have introduced two auxiliary parameters, a G 3? and ip G 9?^. This generic 
equation decomposes the stability problem into two separate parts: one that depends only 
on the individual dynamics and the coupling function, and one that depends only on the 
graph Laplacian and parameter mismatch. Note that the latter not only depends on the 
spectrum of L as in [8], but also on the combination of the eigenvectors and parameter 
mismatch vector. 



Once the stability of Eq. (11 ) is determined as a function of a and the stability of any 



coupled network oscillators as described by Eq. ([TJ, for the given / and H used in Eq. (11), 
can be found by simply setting 

a = gXi (12) 



(13) 



and 

N 

where Xi,Vj^i,6fij can be obtained by the knowledge of the underlying network structure 
L and parameter mismatch pattern. Thus, we have reduced the stability analysis of the 
original miV-dimensional problem to that of an m-dimensional problem with one additional 
parameter, combined with an eigen-problem. 



The associated master stability function (MSF) ^l{a,ilj) of Eq. (11) is defined as: 



lim 



T 



(14) 



when the limit exists, where ^ is a solution of Eq. (11) for the given (a, ijj) pair 



For a given coupled oscillator network by Eq. ([T]), we have the following equation, based 
on the generic MSF Q: 



\ T Jo .^^ 



"1 ^ TV 1 ^ i\ 



rT ^ 



(15) 



j=2 



where Aj are the eigenvalues of the graph Laplacian and ipi is obtained through Eq. (13). 



Thus, once the MSF for the dynamics / and couphng function H has been computed, it 
can be used to compute the asymptotic total distance from single units to the average 
trajectory: < J2iLi ll'W^i(^) ~ ""^jl^)!!^ > [23] for any coupled oscillator network by summing 
up the corresponding Q'^{gXi,ipi) and take the square root. 
In Fig. [T] we plot the MSF for / being Lorenz equations: 

X = a{y — x) 
y = x{r — z) — y 

z = xy — (3z (16) 



[w 



F, z\ 



as individual dynamics in Eq. ([!]) 
10,/? 

H is taken as: H{w) = w, i.e., an identity matrix operator. 



The parameters are chosen as: a 



|, and r is allowed to be adjustable, i.e., r is the /i in Eq. (fl|. The coupling function 



D. Conditions for Stable Near Synchronization 

For near synchronization to appear in the presence of parameter mismatch, it is required 
that the system described by Eq. ([!]) in the absence of parameter mismatch undergoes 
stable identical synchronization, which can be checked by using MSF[8J. In this case, the 
largest Lyapunov exponent of the synchronous trajectory associated with the homogeneous 
variational equation: 

i= [D^f-a-DH]^ (17) 

is negative, and its solution can be written as ^*{t) = $(t, 0)^(0) where $(t, r) is the funda- 
mental transition matrix [21], satisfying 

\mt,T)\\<je-'^'-^^ (18) 

for t > T and some finite positive constants 7 and A. 
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FIG. 1: (Color online) MSF for Lorenz system Eq. (16), with cr = 10, /? = I, adjustable parameter 



being r, and coupling function H being an identity matrix operator. The domain shown here is 
for Q from 5 to 50 and ip from —0.04 to 0.04, while the actual valid domain of MSF could be as 
large as the stability region of the identical Lorenz system (in this case is the region a > Aq where 
Ao is the largest lyapunov exponent of the original Lorenz system (~ 1)). 



Note that the transition to loss of stability at certain time instances can occur due to 
the embedded periodic orbitsjTl [T7], in which case the above inequality will not hold. In 



this paper we consider the situation where Eq. (18) holds for most of the time, with A being 



the Lyapunov exponent of the trajectory associated with Eq. (17), although at certain time 



instances Eq. (18) need not hold, as discussed in [Hill], referred to as bubbling transition [6]. 
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The solution to Eq. (11) can be expressed as: 

e(t) = $(t,0)e(0)+ f ^t,T)biT)dT, 



(19) 



where we have defined 6(r) = D ^^J {s{t) , ii) ■ ip. Under the condition of Eq. (18), we have 
the following bound for ^{t): 



wmw < \m. 



\\m\\+ / \m,r)\\dr.,nv\\h{t)\\ 
Jo t 



< 76 
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7, 



-'*iie(o)ii + ^(i-e 



sup||6(t)|| 



sup ||&(t)|| as t — s> oo. 
A t 



(20) 



Thus, the conditions for stable near synchronization of Eq. ([T]) are: 

1. The corresponding identical system (without parameter mismatch) is stably synchro- 



nized, or equivalently, the associated variational equation Eq. (17) is exponentially 
stable; 



2. The inhomogeneous part 6(r) = D^/(s(r),/i) ■ ip in Eq. (11) is bounded. 



These conditions are sufficient to guarantee the boundness of pairwise distance between any 
two units, so that near synchronous state is stable. 



Eq. (18) and Eq. (19) also allow us to analyze quantitatively the magnitude of asymptotic 



error of a near-identical system such as Eq. ([T|. For all other variables being the same, if the 
magnitude of parameter mismatch is scaled by a factor k, then the corresponding variation 
will become: 



(21) 

where ^ (t) denotes the variation of the original unsealed near-identical system, which follows 



Eq. (19). The first term of both Eq. (19) and Eq. (21) goes to zero according to Eq. (18) 



so that asymptotically the following holds: ^(t) = k ■ C,{t), i.e., the variation is scaled by the 
same factor correspondingly. 

III. EXAMPLES OF APPLICATION 



A. Methodology 



When the units coupling through the network are known exactly, meaning that the pa- 
rameter of each unit is known, then from Eq. (12) and Eq. (13) we may use the fl obtained 
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from MSF at the corresponding (a, V") pairs. In Sec. 3.2 and Sec. 3.3 we illustrate this with 
examples of small networks. 

On the other hand, for large networks, in the case that parameters of individual units 



are not known exactly, but follow a Gaussian distribution: 6fii ~ N{fi,e'^), then in Eq. (10) 
we have: 

TV N 

~ Nm,e') (22) 
assuming the are identical and independent. The standard deviation e may be used, as 



an expected bound for ip in Eq. (13), to compute an expected MSF to predict the possible 
variation of individual units to the average trajectory. In Sec. 3.3 a scale-free network with 
= 500 vertices is used to illustrate. 



In all the examples, the individual dynamics is the Lorenz equation Eq. (16), with param 



eters a = 10,/5 = |, and = 28 + 6ri where 6ri is the parameter mismatch on unit i. The 
coupling function is chosen as H{w) = w with coupling strength g specified differently in 
each example. The variation of individual units to the average trajectory < J2iLi ||'7i(^)|P > 
is approximated by T = 200 with equally time spacing r = 0.01. 



B. Example: Ring Graph 

We consider a small and simple graph to illustrate. The graph as well as three different 
patterns of parameter mismatch are shown in Fig. |2j In Fig. [3] we show the actual variation 
on individual units and that by MSF. 

The MSF predicts well the actual variations found in this near-identical oscillator network, 
in all three cases. Furthermore, the way parameter mismatch are distributed in the graph is 
relevant, as a consequence of Eq. (10). From left to right in Fig. (|2]), the parameter mismatch 



is distributed more heterogeneously, resulting in larger variation along the near synchronous 
trajectory. 
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+ - + 



FIG. 2: (Color online) Ring graph (red circles represent vertices and black lines represent edges) 
with specific parameter mismatch on each unit. The magnitude of parameter mismatch on each 
unit is assumed to be the same, e. The plus/minus sign on a vertex represents the corresponding 
sign of mismatch on that unit, " + " for +e and " — " for — e. So the left graph has the parameter 
mismatch pattern (starting from the top unit): [— e, +e, — e, +e, — e, +e], the middle graph has the 
pattern [— e, — e, +e, — e, +e, +e], and the right graph has the pattern [— e, — e, — e, +e, +e, +e]. 

C. Example: Sequence Networks 

Sequence networks [21] are a special class of networks that can be encoded by the so 
called creation sequence. In Fig. |4]three different sequence networks of the creation sequence 
{A, A, A, B, B, B) under different connection rules are shown. Interestingly, despite the fact 
that the structure of these networks are different, the variation of individual units to the 
average trajectory are the same, under the mismatch pattern [— e, — e, — e, +e, +e, +e], see 
Fig.|5j 

Study on the eigenvector structure on these networks shows that this comes from the 
fact that the eigenvectors of all these three networks are the same, and more importantly, 
the parameter mismatch vector [— e, — e, — e, +e, +e, +e] is parallel to one of the eigenvectors, 
corresponding to the same eigenvalue A = 6 in all three cases. Thus, the only active error 
mode in the eigenvector basis are the same for all three networks, resulting in the same 
variations. 

D. Example: Scale-free Networks 

The synchronization stability of a large network, with the knowledge of the probability 
distribution of parameters, is another interesting problem. To show how an expected MSF 
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FIG. 3: (Color online) Validating MSF on a ring graph. Here the coupling strength is g = 5. The 
units are coupled through a ring graph, with specific parameter mismatch patterns as shown in 
Fig. [2j The vertical axis represents the average variation at each given e. Blue squares, crosses, 
and circles are obtained from actual time series, computed through < J2iLi > where rii{t) 

is the distance from unit i to the average trajectory at time t. Black lines (dashed, solid, and 



dotted) are theoretical prediction \/J2i^i^'^i'^iT'^i) from MSF at {aijipi) paris, where {ai,^|J) are 



computed according to Eq. (12) and Eq. (13). 



will apply, we use a scale-free network as an example. The network is generated using the BA 
model^ : start with a small initial network, consecutively add new vertices into the current 
network; when a new vertex is introduced, it connects to m preexisting vertices, based on 
the preferential attachment rule|9]. The network generated through process is known as a 
BA network, which is one example of a scale-free network. Here we use generate such a BA 
network with = 500 vertices and m = 12. 

In Fig. |6]we show how parameter mismatch affect synchronization on a BA network. The 
parameters on each unit are assumed to follow the Gaussian distribution with mean 28 and 
standard deviation e for each given e. The expected MSF, as described in Sec. 3.1, predicts 
well the actual variation to the average trajectory, see Fig. |6| 
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FIG. 4: (Color online) 2-letter sequence nets[2T] consisting of 6 vertices and 2 layers (red and 
blue circles represent vertices and black lines represent edges) obtained from the same creation 
sequence {A, A, A, B, B, B) under three different rules, on the left the connection rule is ^ ^, 
meaning that whenever a vertex of type B is added into the current net, it connects to all previous 
vertices of type A^ thus a bipartite complete graph is created based on this sequence; on the 
middle the connection rule is i? ^ ^, i3, resulting in a threshold graph; while on the right the rule 
A ^ A^ B; B ^ A, B \s applied to yield a complete graph. Ovals and boxes are used to highlight 
the layer structure: vertices within an oval do not have connections, while vertices within a box 
connect to each other; an thick edge goes from one group to the other connects every vertex in one 
group to all the vertices in the other group. The parameter mismatch pattern here is prescribed 
to coincide with the type of vertices, which is, for the given e: [— e, — e, — e, +e, +e, +e]. 

IV. SUMMARY 

In this paper we analyze the synchronization stability for coupled near-identical oscil- 
lator networks such as Eq. [TJ We show that the master stability equations and functions 
can be extended to this general case as to analyze the synchronization stability. The varia- 
tional equations in the near-identical oscillator case highlight the relevance of eigenvectors as 
well as eigenvalues on the effect of parameter mismatch, which indicates the importance of 
knowledge of the detailed network structure in designing dynamical systems that are more 
reliable. 
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FIG. 5: (Color online) Validating MSF on sequence networks. Here the coupling strength is 5 = 2. 
The parameter mismatch pattern is shown in Fig. |4| The vertical axis represents the average 
variation at each given e. (Blue) markers represent variations obtained from actual time series, 
and (black) dashed lines is the prediction obtained by MSF. Here the MSF line for all three networks 
are the same. 
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